Amplification of Galactic Magnetic Fields 
by the Cosmic-Ray Driven Dynamo 
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ABSTRACT 

We present the first numerical model of the magnetohydrodynamical cosmic-ray (CR) driven 
dynamo of the type proposed by Parker (1992). The driving force of the amplification pro- 
cess comes from CRs injected into the galactic disk in randomly distributed spherical regions 
representing supernova remnants. The underlying disk is differentially rotating. An explicit re- 
sistivity is responsible for the dissipation of the small-scale magnetic field component. We obtain 
amplification of the large-scale magnetic on a timescale 250 Myr. 

Subject headings: Galaxies: ISM, Magnetic Fields; ISM: Cosmic Rays, Magnetic Fields; MHD: Dynamos 



1. Introduction 

In 1992 Parker discussed the possibility of a new 
kind of galactic dynamo driven by galactic CRs 
accelerated in supernova remnants. This dynamo 
contains a network of interacting forces: the buoy- 
ancy force of CRs, the Coriolis force, the differen- 
tial rotation and magnetic reconnection. Parker 
estimated that such a dynamo is able to amplify 
the large scale magnetic field on timescales of the 
order of 10 8 yr. 

Over the last decade we have investigated the 
different physical properties and consequences of 
Parker's idea and scenario by means of analyti- 
cal calculations and numerical simulations (see eg. 
Hanasz & Lesch 1998, Lesch & Hanasz 2003 and 
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references therein) . Here we present the first com- 
plete magnetohydrodynamical three-dimensional 
simulation including the full network of relevant 
interacting mechanisms. 

It is the aim of our contribution to show that 
Parker's CR driven dynamo indeed acts efficiently 
on timescales comparable with the disk rotation 
time. In the next two Sections we describe the 
physical elements of the model and the system of 
equations used in numerical simulations. Section 4 
presents the numerical setup, Sections 5 and 6 in- 
form the reader about the results on the structure 
of the interstellar medium including CRs and mag- 
netic fields, the strength of the amplified magnetic 
field and the spatial structure of the mean mag- 
netic field. We summarize our results very briefly 
in Section 7. 

2. Elements of the model 

We performed computations with the aid of the 
Zeus-3D MHD code (Stone & Norman 1992a,b), 
which we extended with the following features: 

(1) The CR component, a relativistic gas de- 
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scribed by the diffusion-advection transport equa- 
tion (see Hanasz & Lesch 2003b for the details 
of numerical algorithm). Following Jokipii (1999) 
we presume that CRs diffuse anisotropically along 
magnetic field lines. (2) Localized sources of CRs: 
supernova remnants, exploding randomly in the 
disk volume (see Hanasz & Lesch 2000). (3) Resis- 
tivity of the ISM (see Hanasz et al. 2002, Hanasz 
& Lesch 2003a) responsible for the onset of fast 
magnetic reconnection (in this paper we apply the 
uniform resistivity) . (4) Shearing boundary condi- 
tions and tidal forces, following the prescription by 
Hawlcy, Gammie & Balbus (1995), aimed to model 
differentially rotating disks in the local approxima- 
tion. (5) Realistic vertical disk gravity following 
the model of ISM in the Milky Way by Ferriere 
(1998). 

3. The system of equations 

We apply the following set of resistive MHD 
equations 

^ + V-( P V) = 0, (1) 
g + V.(eV) = -p(V-V), (2) 

_ + (V.V)V = --v( p + Pct + -) 
B * VB 

+— A 2fl x v + 2qn 2 xe^ + g z e z , (3) 

47T/9 

— = Vx(VxB) + V AB, (4) 

p=(7-l)e, 7-5/3 (5) 

where q — — d lnil/d InR is the shearing parame- 
ter, (R is the distance to galactic center), g z is the 
vertical gravitational acceleration, -q is the resistiv- 
ity, 7 is the adiabatic index of thermal gas, the gra- 
dient of CR pressure Vp cr is included in the equa- 
tion of motion (see Hanasz & Lesch 2003b) and 
other symbols have their usual meaning. The uni- 
form resistivity is included only in the induction 
equation (see Hanasz et al. 2002). The adopted 
value J] — 1 exceeds the numerical resistivity for 
the grid resolution defined in the next section (see 
Kowal et al. 2003). The thermal gas component 
is currently treated as an adiabatic medium. 



The transport of the CR component is de- 
scribed by the diffusion-advection equation 

(e cr V) = V (k Ve cr ) - Pcr (V • V)+Q SN , 

(6) 

where Qsn represents the source term for the CR 
energy density: the rate of production of CRs in- 
jected locally in SN remnants and 

Per = (7cr - l)e cr , 7cr = 14/9. (7) 

The adiabatic index of the CR gas 7 cr and the 
formula for diffusion tensor 

Kij = K±Sij + (K\\-K±)ninj, rii = Bi/B, (8) 

are adopted following the argumentation by Ryu 
et al. (2003). 

4. Numerical simulations 

We performed numerical simulations in a 3D 
Cartesian domain 500pc x lOOOpc x 1200pc , 
extending symmetrically around the galactic mid- 
plane from z — — 600pc up to z — 600pc , with 
the resolution of 50 x 100 x 120 grid zones in di- 
rections x, y and z, corresponding locally to cylin- 
drical coordinates r, (f> and z, respectively. The 
applied boundary conditions are periodic in the y- 
direction, sheared-periodic in the x direction and 
outflow in the z direction. The computational vol- 
ume represents a 3D region of disk of a galaxy 
similar to the Milky Way. 

The assumed disk rotation is represented locally 
by the angular velocity ft = 0.05Myr _1 and by a 
flat rotation curve corresponding to q = 1. We ap- 
ply the vertical gravity profile determined for the 
Solar neighborhood (see Ferriere 1998 for the for- 
mula). We assume that supernovae explode with 
the frequency 2kpc _2 Myr _1 , and assume that 10 
% of the 10 51 erg kinetic energy output from SN 
is converted into the CR energy. The CR en- 
ergy is injected instantaneously into the ISM with 
a Gaussian radial profile (r SN — 50pc) around 
the explosion center. The explosion centers are 
located randomly with a uniform distribution in 
the x and y directions and with a Gaussian dis- 
tribution (scaleheight H — lOOpc) in the vertical 
direction. The applied value of the CR parallel 
diffusion coefficient is K\\ = 10 4 pc 2 Myr^ 1 = 3 x 
10 27 cm 2 s _1 (i.e. 10 % of the realistic value) and 
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Fig. 1. — Slices through the computational volume at t = 2000Myr . Panels (a),(b) and (c) show CR energy 
density with vectors of magnetic field in yz, xz and xy planes respectively, panels (d), (e) and (f) show gas 
density with velocity vectors in the same planes. 
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Fig. 2. — Evolution of magnetic field in the computational volume. The left panel shows the total mag- 
netic energy normalized to the initial gas thermal energy (dash-dot line) along with normalized energies of 
azimuthal magnetic field (full line), radial magnetic field (dotted line) and vertical magnetic field (dashed 
line). The right panel shows the evolution of azimuthal magnetic flux x(— 1) (full line) and radial magnetic 
flux xlO (dotted line), both normalized to the initial value of the azimuthal magnetic flux. 
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Fig. 3. — Vertical structure of the mean magnetic field for a sequence of time instants. Spatial averaging is 
done over planes z=const. Full line represents the mean azimuthal magnetic field and dotted line the radial 
magnetic field multiplied by 10. The amplification effect is reflected in the varying range of the horizontal 
axis. 



the perpendicular one is Kj_ = 10 3 pc 2 Myr 1 = 
3 x 10 26 cm 2 s" 1 . 

The initial state of the system is a magnetohy- 
drostatic equilibrium with a horizontal purely az- 
imuthal magnetic field of the strength correspond- 
ing to Pmag/pgas = 10~ 8 . The initial CR pressure 
in the initial state is equal to zero. The initial 
gas density at the galactic midplanc is 3 H atoms 
cm~ 3 and the initial isothermal sound speed is 
c s ; = 7km s _1 . 

5. Structure of interstellar medium result- 
ing from the CR-MHD simulations 

In Fig. 1 we show the distribution of CR gas 
together with magnetic field, and thermal gas den- 
sity together with and gas velocity in the compu- 
tational volume at t = 2000 Myr. 

One can notice in panel (a) a dominating hor- 
izontal alignment of magnetic vectors. The CR 
energy density is well smoothed by the diffusive 
transport in the computational volume. The ver- 
tical gradient of the CR energy density is main- 
tained by the supply of CRs around the equatorial 
plane in the disk in the presence of vertical grav- 
ity. In panel (c) one can notice that at the height 
z = 370pc the dominating magnetic vectors are 
inclined with respect to the azimuthal direction, 
i.e. the radial magnetic field component is on av- 
erage about 10 % of the azimuthal one. 

The CR energy density is displayed in units in 
which the thermal gas energy density correspond- 



ing to p = 1 and the sound speed c s i = 7km s _1 is 
equal to 1. We note that the CR energy density 
does not drop to zero at the lower and upper z 
boundaries due to our choice of outflow boundary 
conditions for the CR component. We note also 
that almost constant mean vertical gradient of CR 
energy density is maintained during the whole sim- 
ulation. 

The velocity field together with the distribution 
of gas density, is shown in the next panels (d), (e) 
and (f). The shearing pattern of velocity can be 
noticed in the horizontal slice (f). The vertical 
slices (d) and (e) show the stratification of gas 
by the vertical gravity, acting against the vertical 
gradients of thermal, CR and magnetic pressures. 

6. Amplification and structure of the 
mean magnetic field 

In the following Fig. 2 we show how efficient 
is the amplification of mean magnetic field re- 
sulting from the continuous supply of CRs in su- 
pernova remnants. First we note the growth of 
the total magnetic energy, by 7 orders of magni- 
tudes during the period of 2 Gyr. Starting from 
t ~ 300 Myr the growth of magnetic energy repre- 
sents a straight line on a logarithmic plot, which 
means that the magnetic energy grows exponen- 
tially. The c-folding time of magnetic energy de- 
termined for the period t = 400 -j- 1500 Myr is 115 
Myr. Around t = 1500 Myr, the growth starts to 
slow down as the magnetic energy approaches an 
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equipartition with the gas energy. 

The other three curves in the left panel of Fig. 2 
show the growth of energy of each magnetic field 
component. It is apparent that the energy of ra- 
dial magnetic field component is almost an order 
of magnitude smaller than the energy of vertical 
magnetic field component which is almost one or- 
der of magnitude smaller than the energy of the 
azimuthal one. This indicates that the dynam- 
ics of the system is dominated by the buoyancy 
of CRs and that magnetic reconnection efficiently 
cancels the excess of the random magnetic fields. 

In the right panel of Fig. 2 we show the time 
evolution of the normalized, mean magnetic fluxes 
®x(t)/® y (t = 0) and $ y (t)/®y(t = 0), where 
<& x (t) and & y (t) are respectively magnetic fluxes at 
moment t, threading vertical planes perpendicular 
to x and y axes respectively, and averaging is done 
over all possible planes of a given type. We find 
that the radial magnetic flux <& x starts to deviate 
from zero, as a result of Coriolis force and open 
boundary conditions. Due to the presence of dif- 
ferential rotation the azimuthal magnetic field is 
generated from the radial one. The azimuthal flux 
grows up by a factor of 10 in the first 800 Myr of 
the system evolution and then drops suddenly, re- 
verts and continues to grow with the opposite sign 
undergoing amplification by more three orders of 
magnitudes, with respect to the initial value. 

In order to examine the structure of the mean 
magnetic field we average of B x and and B y across 
constant z-planes. The results are presented in 
Fig. 3 for t = (the initial magnetic field) and 
then for t = 500, 1000, 1500, 2000 and 2300 Myr. 
We find that the mean magnetic field grows by a 
factor of 10 within about 500 Myr, which gives an 
e-folding time close to 250 Myr. We note that an 
apparent wavelike vertical structures in (B x ) and 
(By) formed from the initial purely azimuthal, uni- 
directional state of B y and B x = 0. The evolved 
mean magnetic field configuration reaches a qua- 
sisteady pattern which is growing in magnitude 
with apparent vertical reversals of both compo- 
nents of the mean magnetic field. We note also 
that the magnetic field at the disk midplane re- 
mains relatively weak. 

A striking property of the mean magnetic field 
configuration is the almost ideal coincidence of 
peaks of the oppositely directed radial and az- 
imuthal field components. This feature corre- 



sponds to a picture of an a — 51-dynamo: the az- 
imuthal mean magnetic component is generated 
from the radial one and vice versa. 

In order to understand better what kind of dy- 
namo operates in our model, we computed the y- 
component of the electromotive force (£mf, y ) = 
(v z B x — v x B z ), averaged over constant z-plancs 
and checked that d(B x )/dt ~ —d(£ m ^ y )/dx with 
a reasonable accuracy. However, we found that 
the space averaged £ m f, y fluctuates rapidly in time, 
so that the approximation of (£ m f, v ) by a yy (B y ) 1 
(where a yy is a component of the fluid helicity ten- 
sor), implies that a yy oscillates rapidly in time. 
This property points our model toward the inco- 
herent a — il dynamo described by Vishniac & 
Brandenburg (1997). Finally, we checked that the 
magneto-rotational instability (Balbus and Haw- 
ley 1991) does not seem to play a significant role 
in our dynamo model. Due to the weakness of 
the initial magnetic field, the wavelength of most 
unstable mode of this instability remains shorter 
than the cell size for the first half of the simulation 
time. 

7. Conclusions 

We have described the first numerical experi- 
ment in which the effect of amplification of the 
large scale galactic magnetic field was achieved 
by the (1) continuous (although intermittent in 
space and time) supply of CRs into the interstel- 
lar medium, (2) shearing motions due to differen- 
tial rotation and (3) the presence of an explicit 
resistivity of the medium. 

We observed in our experiment the growth of 
magnetic energy by seven orders of magnitude and 
the growth of magnetic flux by a factor of 1300 in 
2150 Myr of the system evolution. We found that 
the large scale magnetic field grows on a timescale 
250 Myr, which is close to the period of galactic 
rotation. 

Therefore the galactic dynamo driven by CRs 
appears to work very efficiently, as it was sug- 
gested by Parker (1992). It is a matter of future 
work to verify whether the presented model is a 
fast dynamo, i.e. whether it works with a similar 
efficiency in the limit of vanishing resistivity. 

This work was supported by the Polish Com- 
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